% This file analyzes the experiments intended to reveal the tensor monopole
% generates Fig2A

alpha_fit=linspace(0,pi/2,41);
load('./data/Fig2A.mat')

% g components
cc =[0.0549    0.2706    0.6275
    0.9020         0    0.0627
    0.0980    0.5412    0.1843
    0.8941    0.2314    0.5333
    0.9647    0.3294    0.0510
    0.0941    0.5294    0.9098
    0.5686    0.7490    0.1020
    0.9882    0.7098    0.0588
    0.1020    0.5686    0.4745];
figure
hold on
for hlp=1:6
    switch hlp
        case 1
            plot(alpha_fit/pi*180,ones(size(alpha_fit))/2,'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off')
        case 2
            plot(alpha_fit/pi*180,1/4*cos(alpha_fit).^2.*(2-cos(alpha_fit).^2),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off')
        case 3
            plot(alpha_fit/pi*180,1/4*sin(alpha_fit).^2.*(2-sin(alpha_fit).^2),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off')
        case 4
            plot(alpha_fit/pi*180,zeros(size(alpha_fit))/2,'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off')
        case 5
            plot(alpha_fit/pi*180,zeros(size(alpha_fit))/2,'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off')
        case 6
            plot(alpha_fit/pi*180,1/32*(cos(4*alpha_fit)-1),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off')
    end
end
for hlp=1:6
    errorbar(alpha_sweep/pi*180,g_munu_save(hlp,:),g_munu_err_save(hlp,:),'d','MarkerSize',12,'LineWidth',2,'HandleVisibility','off','Color',cc(hlp,:),'HandleVisibility','on')
end
hold off
xlim([0,90]);set(gca,'XTick',[0:30:90])
xlabel('\alpha ({\circ})','FontSize',18)
ylabel('g_{\mu\nu}','FontSize',18)
set(gca,'FontSize',18)
lgd=legend('g_{\alpha\alpha}','g_{\beta\beta}','g_{\phi\phi}','g_{\alpha\beta}','g_{\alpha\phi}','g_{\beta\phi}');lgd.FontSize=18;
set(lgd,...
    'Position',[0.1452 0.5873 0.7446 0.0798],...
    'Orientation','horizontal',...
    'FontSize',18);
% set(lgd,'Location','best','Orientation','horizontal','FontSize',18);

function [rabi_g,E_total,full_map,g_munu_save] = QGT_rabi_freq_detune(alpha_sweep,beta,phi,rd,r,flag)

rabi_g = zeros(2,length(flag),length(alpha_sweep));
E_total = zeros(2,length(alpha_sweep)); % 1: SQ, 2: DQ
full_map = zeros(1,length(alpha_sweep));
dpert=0.005;
g_out = zeros(length(flag),length(alpha_sweep));
g_munu_save = zeros(6,length(alpha_sweep));

for hlp=1:length(alpha_sweep)
    theta=alpha_sweep(hlp);
    H0=Ham_detune(theta,beta,phi,rd);
    [V,E]=eig(H0);
    [E,idx]=sort(diag(E));
    
    E_total(:,hlp)=[abs(E(2)-E(1))*2;abs(E(3)-E(1))]*r; % note we already take 2X the SQ frequency here
    
    V=V(:,idx);
    um=V(:,1)/sqrt(V(:,1)'*V(:,1));
    u0=V(:,2)/sqrt(V(:,2)'*V(:,2));
    up=V(:,3)/sqrt(V(:,3)'*V(:,3));
    
    for hlp_flg=1:length(flag)
        flg_tmp=flag{hlp_flg};
        
        theta_oscillate=flg_tmp(1);
        beta_oscillate=flg_tmp(2);
        phi_oscillate=flg_tmp(3);
        
        for hlp_omega=1:2 % omega
            dH=Ham_detune(theta+dpert*theta_oscillate,beta+dpert*beta_oscillate,phi+dpert*phi_oscillate,rd);
            dH=(dH-H0)/dpert;
            
            switch hlp_omega
                case 1
                    rabi_g(hlp_omega,hlp_flg,hlp)=abs(um'*dH*u0*r);
                    
                    g_out(hlp_flg,hlp)=g_out(hlp_flg,hlp) + ...
                        abs(um'*dH*u0*r)^2/((E(1)-E(2))*r)^2;
                    
                case 2
                    rabi_g(hlp_omega,hlp_flg,hlp)=abs(um'*dH*up*r);
                    
                    g_out(hlp_flg,hlp)=g_out(hlp_flg,hlp) + ...
                        (rabi_g(hlp_omega,hlp_flg,hlp))^2/((E(1)-E(3))*r)^2;
            end
        end
    end
    
    g_munu_save(1,hlp)=g_out(1,hlp);
    g_munu_save(2,hlp)=g_out(2,hlp);
    g_munu_save(3,hlp)=g_out(3,hlp);
    g_munu_save(4,hlp)=(g_out(4,hlp)-g_out(5,hlp))/4;
    g_munu_save(5,hlp)=(g_out(6,hlp)-g_out(7,hlp))/4;
    g_munu_save(6,hlp)=(g_out(8,hlp)-g_out(9,hlp))/4;
    
end
end

function Ham_out = Ham_detune(alpha,beta,phi,rd) %
Ham_out=[0,cos(alpha)*exp(-1i*beta),0;cos(alpha)*exp(1i*beta),0,sin(alpha)*exp(1i*phi);...
    0,sin(alpha)*exp(-1i*phi),0]+diag([rd,0,-rd]);
end

function [QT,QT_err] = Herror_charge(alpha_sweep,H_berry_curvature, H_berry_err)
% calculates the topological charge
% takes alpha_sweep which is not equidistant

QT=sum((H_berry_curvature(1:end-1)+H_berry_curvature(2:end)).*(alpha_sweep(2:end)-alpha_sweep(1:end-1)));

QT_err=sqrt(sum((H_berry_err(1:end-1).^2+H_berry_err(2:end).^2).*(alpha_sweep(2:end)-alpha_sweep(1:end-1)).^2));
end

